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PACS 82 . 70 . Dd - Colloids 

PACS 64 . 70 . D- - Solid-liquid transitions 

PACS 61 . 20 . Ja - Computer simulation of liquid structure 

Abstract. - Using grand-canonical Monte Carlo simulations we investigate the impact of charged 
walls on the crystallization properties of charged colloidal suspensions confined between these 
walls. The investigations are based on an effective model focussing on the colloids alone. Our 
results demonstrate that the fluid-wall interaction stemming from charged walls has a crucial 
impact on the fluid's high-density behavior as compared to the case of uncharged walls. In 
particular, based on an analysis of in-plane bond order parameters we find surface-charge-induced 
freezing and melting transitions. 



' ' Introduction. — The ordering of charged nanoparti- 
y cles (colloids) at charged surfaces and in confining geome- 
q tries is a topic receiving continuous and intense attention 
O since decades. Such systems occur in a wealth of con- 
~ taxts such as in material science (e.g., as candidates for 

fNiijovel self- assembled structures [TJH]), in microfluidics in- 
£>volving charged channels [3J, in biological contexts (e.g., 

OC proteins at charged lipid membranes [HE]) and in the 
"^ physics of dusty plasmas [6j-[8]. Moreover, charged colloids 
Ijave turned out to be excellent model systems with tun- 



in 



able (screened electrostatic) interactions, which are ideally 
suited to study confinement-related equilibrium and non- 



f-*. equilibrium phenomena such as freezing in unusual crystal 
^H Structures [8HT4] , oscillatory depletion forces [15l[T6] , and 
!* shear- induced distortions of colloidal films [T7J[T8] ■ 

JT Despite intense research both from the experimental and 
5-j from the theoretical (statistical-mechanical) side, we are 
still far away from a full understanding of the interplay 
of macroion charges, counterions, salt and surface charges 
and their precise role for the physics of confined charged 
suspensions. One prominent, long-standing problem is the 
occurrence of attractive forces between like-charged col- 
loids between glass surfaces [T5H55] , a phenomenon which 
has been attributed to strong countcrion correlations (typ- 
ical for multivalent counterions) or countcrion depletion 
|23j . Another topic, which is also relevant for moderately 
coupled systems, concerns the impact of wall charges on 
the microscopic concentration profiles [55] ■ Indeed, recent 
experimental studies [1111571155] based on x-ray diffraction 
methods have reported novel effects such as salt-induced 



trapping of colloidal monolayers [57] and surface-charge- 
induccd stabilization of layers with respect to buckling 
transitions [28] . 

In the present letter we report computer simulation re- 
sults on yet another phenomenon, that is, the role of sur- 
face charges on the lateral ordering of a charged colloidal 
suspension between two planar, like-charged walls. In an 
earlier, combined theoretical-experimental study |26] we 
have already shown that such wall charges (which stem 
from the countcrion dissociation into the adjacent fluid) 
have a crucial impact on the effective fluid-solid interac- 
tion and the resulting oscillatory structural forces. More- 
over, the theoretically predicted effects for low densities 
coincide with experimental results |26| . 

Here we go one step further and investigate the confined 
fluid's crystallization properties. As in our previous stud- 
ies [m[roi2"rjll2"M!n"j we employ grand canonical Monte 
Carlo (GCMC) simulations based on an effective model 
of the confined colloidal suspension. Our present results 
demonstrate, for the first time, that surface charges do not 
only affect the vertical ordering [55], but also the (cap- 
illary) freezing transitions of the system. In particular, 
variation of the surface charge can induce freezing and 
melting transitions, and it can even change the character 
of the solid phases. 

Model system. — Following our previous work [TT1 
[T51[5rjl[5M5I] we use a model based on Derjaguin-Landau- 
Verwey-Overbeek (DLVO) theory [35], where the inter- 
action between the charged colloidal spheres suspended 
in a solution is described on an effective level assuming 
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relatively small macroion density and charge. All other 
contributions from the dissociated (monovalent) counteri- 
ons and the salt ions are treated implicitly. The resulting 
interaction potential reads 



wdlvoM = W 



exp(— nr) 



with the prcfactor 



W 



{Ze f 

Atteqe 



cxp(ftcr) 



(1) 



(2) 



In eq. (J2J, o is the diameter of the particles, eo is the 
elementary charge, eo and e are the permittivities of the 
vacuum and the solvent, respectively, and Z = Z/(l + 
na/2) is an effective valency assuming that the macroion 
density is small compared to the salt concentration. The 
valency of the particles is Z . The inverse Debye screening 
length k appearing in cq. ([T]) is defined as 



e ek B T 



(Zp+2IN A ), 



(3) 



with macroion density p, ionic strength /, Avogadro's con- 
stant N A , Boltzmann constant k B , and temperature T. 
Thus, we do not take into account any effect of the wall 
counterions on the effective fluid-fluid interaction poten- 
tial. A justification is given in the discussion of fig. [3] 

As in our previous work [TTI1 6, 26, 29-151] we supplement 
the DLVO interaction by a soft-sphere potential ugg(r) = 
4e ss (a/r)_ 12 . 

The fluid-solid interaction between macroions and walls 
based on linearized Poisson-Boltzmann (PB) theory for 
small Debye screening lengths and the linear superposition 
approximation (LSA) reads [3U1IMJ 
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k B T 
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exp 



(-K-w(z)(z- -)j J 



(4) 

where 7f/s = tanh (eoipp/s/4:k B T), and V'f/s is the sur- 
face potential of the fluid particles (F) and the solid walls 
(S). Contrary to |34| . however, the screening parameter 
Kyf(z) depends in our model on the ^-coordinate, reflect- 
ing the inhomogeneous distribution of counterions released 
by the walls. Specifically, using the exact PB expression 
for the counterion profile [23] in the dilute limit, and per- 
forming an averaging procedure, our expression for the 
space-dependent wall screening parameter reads |26] 



Kw(z) 



e ek B T 



Zp- 



21 N A + ^ 
e z 



(5) 



Fig. 1: (a) Fluid-solid interaction potential Up S (z) and (b) local 
density profiles p*(z) of the macroion distribution for L z = 
2.4ct and p£ = 0.61 at surface potentials |^>g| = (solid lines), 
|V>s| = SOmV (dashed lines), and |^ig| = 160?nV (dot-dashed 
lines). The confining walls are located at z = — 1.2a and 1.2a. 



resulting from a competition between Coulomb repulsion 
and screening [26j . The same qualitative behavior of ufs 
is observed with a homogeneous ansatz for k\v |13U23] , as 
discussed in detail in 0S]- We note, that our ansatz in 
cq. (JS|) yields a dependence of resulting surface forces on 
the wall charge consistent with the experiment. Indeed, 
the data in fig. [1^, arc based on realistic values of ips such 
as silica (|"0s| = 80mV) and mica (IV'sl = 160™^) which 
are typical materials for experiments |26j . For increas- 
ing erg the repulsion range of the surfaces increases and 
the pore is effectively narrowed. Above a certain charge 
the wall screening begins to dominate the Coulomb re- 
pulsion and yields an effective broadening. However, the 
effective pore width for the highest surface charge con- 
sidered (mica) is still smaller compared to the uncharged 
case. The non-monotonic behavior of ufs has a strong 
impact on the particle layering of the macroions as shown 
in fig. [TJ) where we present density profiles p z |11| . We 
note that these effects remain when k (eq. ([3])) appearing 
in the fluid-fluid interaction is replaced by a screening pa- 
rameter involving the surface charge explicitly, as we show 
in fig. [3] Our model only contains pure repulsive interac- 
tions. As motivated in [26] like-charge attraction [19 25 
between the macroions is not expected in the parameter 
range of interest. 

Parameters are set in order to mimic the conditions of 
Colloidal-Probe Atomic- Force-Microscope experiments in- 
volving dilute silica suspensions with low salt concentra- 
tions [2"6"ll2"9"H3"T] . Resulting values of the coupling and 
screening parameters are W* = (k B T)~ 1 W ~ 80 — 90 
and k* = k<t k 2.5 — 2.7 [IT]. In-plane correlation 
functions <7||(r||) = (iVi(j"||, Ar)) / (Nipi2irr\\Ar) for a cer- 
tain layer I with N\ particles within a small interval 
Ar for the in-plane particle separation j-h = \fkj\ = 



\/( x kj +ylj) °f P ar ticles k and j, and the particle area 
density p\ are extracted from the simulations. Another 
key quantity is the local bond order parameter ip n = 

tential'wFsOzj = k B Tu* FS (z) = u^f(z) + u|| A (a;) shown in (V^W* EkLT l / N k Sji cxp(in6» fcj ) j with the angle 
fig.[TJi reveals a non-monotonic behavior as function of o~s 9kj enclosing the bond vector fkj and an arbitrary in-plane 



with the surface charge density erg. The fluid-solid 
interaction is supplemented by a soft-wall repulsion 
i%Y(z) = (4/5)7re w (cr/z) 9 . The resulting fluid-solid po- 
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axis. It is a meassure for the local in-plane order of a cer- 
tain layer parallel to the confining walls with iVi ayer parti- 
cles [TT] . In order to find the N^ next neighbors of a par- 
ticle we employed a Delaunay triangulation |35j . Crystal- 
like ordering is signalled by values of tp n larger than 0.6 as 
motivated by the results in fig. [B] 

Results and Discussion. — We have performed 
GCMC simulations with confined systems containing 
about 500, 1000 and 2000 particles. After equilibration of 
the systems using 50, 000 - 500, 000 MC steps we collected 
averages of relevant physical quantities with 100, 000 - 
1, 000, 000 MC steps. Our present investigations focus on 
nanoconfined suspensions characterized by surface separa- 
tions L z yielding one- or two layer systems. 

A first impression of the role of both, the value of L z 
and the surface charge, is given in fig. [31 where we plotted 
bond order parameters as function of L z for p£ = 0.61. 
Specifically, we consider ?/>6 and ^8 which are sensitive to 
hexagonal and squared order, respectively. We used tps 
instead of ip4 since the next neighbors were determined 
via Delaunay triangulation where on average the particles 
have six next neighbors. Therefore, in contrast to ij)% the 
parameter ^4 cannot become 1.0 even in a perfect squared 
lattice. As known from previous work |11) the bond or- 
der parameters exhibit an alternating behavior as function 
of L z stemming from a competition between layering and 
in-plane ordering. As seen in fig. [5] this behavior is signif- 
icantly affected by the wall potential. On one hand the 
phase of the alternation is shifted forwards larger sepa- 
rations when the uncharged walls are replaced by silica 
walls since the interaction range of the wall potential in- 
creases (see fig. [T]). On the other hand the height of the 
maxima of tpg and "0s becomes smaller because the layers 
are softened (26 . For larger surface charges up to mica 
(IV'sl ^ lOOmV) the opposite behavior occurs. The oscil- 
lations are shifted back and the maxima are enhanced. In 
particular, looking at a certain separation such as, e. g., 
L z = 2.4cr the parameter ipg jumps from 0.5 (uncharged) 
to 0.7 (silica) and back to 0.3 (mica). Similar behavior is 
observed for t/jg. Mica surfaces exhibit a bond order larger 
than 0.6 whereas the other wall types yield much smaller 
values. 

This charge-induced ordering and disordering is also vis- 
ible from the in-plane correlation functions 511 (rii). In 
fig. [3^, correlation functions arc plotted for the three differ- 
ent wall types, L z = 2.4<7 and p^ = 0.61. For uncharged 
surfaces the correlations behave still fluid-like. On the 
other hand, silica and mica walls both induce crystal-like 
correlation functions reflecting hexagonal and squared or- 
der for silica and mica, respectively. Data for L z = 2.9cr 
are shown in fig.[3}D. In this case, uncharged and mica walls 
exhibit crystal-like translational correlations whereas silica 
walls correspond to fluid-like behavior. Both results are 
consistent with the charge-induced ordering of the parti- 
cles seen in fig. [2] 

At this point, we briefly discuss the influence of the ad- 
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Fig. 2: Bond order parameters (a) ipe and (b) ijjs as function of 
the wall separation L z at p£ = 0.61 for uncharged (solid line), 
silica (dashed line), and mica (dot-dashed line) walls. The 
triangles and squares are the simulation data, and the lines are 
a guide to the eye. 



ditional wall counterions on the effective fluid-fluid inter- 
action. Such an approach was suggested in [23j , where the 
screening parameter of the DLVO fluid-fluid interaction 
includes a counterion contribution p cw stemming from the 
surface charge defined as an average of the microscopic 
(PB) profile, that is, 



K = k(zi,z 2 ) = 



^{Zp + 2IN A + p cw ( Zl ,z 2 )). 

e ek B -l- 



(6) 
We tested the importance of this modification of the fluid- 
fluid interaction by calculating corresponding in-plane cor- 
relations at L z = 2.9cr. The results in fig. [3Jd, however, 
show that the effect is marginal. 

We now consider the impact of wall charges on the sys- 
tems ordering behavior as function of the chemical poten- 
tial (and thus, the density), focussing on the wall separa- 
tion L z = 2.4cr. Figures 0^ and @J) show corresponding 
results for tp 6 and -0s for uncharged, silica (\ips\ = 80m V), 
and mica walls (IV'sl = 160mV). Whereas tp6 for un- 
charged and mica walls is less than 0.5, silica walls induce 
a significant hexagonal ordering which becomes particu- 
larly pronounced for p* > 40.5 (p£ > 0.59). Indeed, at 
the largest p* considered the order parameter ipg has val- 
ues of about 0.75 consistent with the clear hexagonal or- 
dering seen in fig. [4J3. Another type of charge-induced 
ordering occurs for mica walls. This is clearly visible 
from fig. [3}3, which reveals a marked increase of ip% for 
mica walls, whereas corresponding order parameters re- 
main constant or even decrease for uncharged and silica 
walls, respectively. Summarizing, for large enough p* an 
increase of the surface potential from |^g| =0 (uncharged) 
to jf/'sl = 80mV (silica) and | Vs I = WOmV^ (mica) yields a 
change of the in-plane order from non-ordered to hexago- 
nal and finally to squared. Moreover, the onset of transla- 
tional ordering is slightly shifted towards smaller chemical 
potentials when silica is replaced by mica (uncharged walls 
do not order at all in the density range considered). 

To confirm the picture that the charged walls lead to 
true crystallization phase transitions we have computed 
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Fig. 3: (a) In-plane correlation functions for L z = 2.4<j using 
the wall types uncharged (inset, thick black solid line), silica 
(thin black solid line) and mica (thin grey solid line), (b) In- 
plane correlation functions for L~ = 2.9a with wall types as in 
(a). Included are results involving the wall-counterion contri- 
bution in the fluid-fluid interaction for silica (thin black dashed 
line) and mica (thin grey dashed line). 



fluctuations (susceptibilities) of the bond-order parame- 
ters defined as Xn = NUipn) ~ (VO )■ These fluctua- 
tions are expected to increase near a phase transition, 
i. e., at the phase transition % n should exhibit a maxi- 
mum [36]. Figure 2fc contains the susceptibility xe corre- 
sponding to the data for ijjq in fig. |3Ji. This quantity pos- 
sesses for silica a maximum around /i* = 41.5 (p£ — 0.61) 
which can be identified with the onset of hexagonal order- 
ing. As expected x§ remains small for mica. Considering 
squared structures xs exhibits a maximum at fj,* = 39.5 
(/?£ = 0.58) in fig. |4ji which is consistent with the strong 
increase of -08 m fig- Et>- Although the bond order pa- 
rameter 08 remains small for silica walls it also shows a 
maximum in xs- This non-expected maximum for silica in 
fig. 0Ji is related to strong bond order fluctuations when 
the hexagonal order in the system emerges (see ipe and 
Xe)- 

We note in this context that the results for the sus- 
ceptibilities close to the phase transition are plagued by 
the fact that the simulated system "jumps" between two 
(coexisting) phases, despite the large free energy barrier 
separating fluid and solid phases. Indeed, most of the 
data points in figs. 0J; and 0J1 do not contain any jumps, 
an exception being the data corresponding to /j,* = 41.52 
and jU* = 42.48 for silica walls with large X6 arid xs- For 
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Fig. 4: Bond angle order parameters (a) i/>6 and (b) Vs with 
corresponding susceptibilities (c) x& an d (d) x& as function of 
^* for uncharged (solid line), silica (dashed line), and mica 
(dot-dashed) walls with a separation of i z = 2.4<r. The trian- 
gles and the squares are simulation results, whereas the lines 
are a guide to the eye. Corresponding top-view snapshots for 
jU* = 41.52 and (e) silica walls (hexagonal phase) and (f) mica 
walls (squared phase). 



/i* = 41.52 we show in fig. [5^ results for the order param- 
eter distribution P(ipe). The latter reflects a double-peak 
structure typical for a first-order transition near coexis- 
tence [SUES]- The peaks overlap due to the similarity of 
the order parameter in the two phases. The corresponding 
behavior of the instantaneous value of tpg as function of 
MC "time" clearly shows not only fluctuations on short 
scale, but also jumps between the two phases on a longer 
time scale. 

The first-order character of the observed phase transi- 
tions is also indicated by the behavior of the susceptibil- 
ities as functions of the system size. Specifically, for a 
second-order phase transition the susceptibility at phase 
coexistence is supposed to diverge in the thermodynamic 
limit whereas for a first-order transition this maximum 
docs not depend on N [36] . As an example, we consider 
in fig. [5Jd results for the emerging squared phase between 
mica surfaces (see figs. Up and Hpl). Computing xs f° r 
the particle numbers N = 500, 1000, 2000 we see that the 
maximum remains more or less the same (fig. \5jp). The 
small deviations are supposed to be due to statistical er- 
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Fig. 5: (a) Distribution P(ipe) of the bond angle order tpe for 
L z = 2.4<r and p£ = 0.61 (/x* = 41.52) using silica walls. The 
inset shows the actual non-averaged values of ipg during the 
MC simulation in the averaging interval, (b) System size de- 
pendence of the bond order fluctuations \s as function of the 
chemical potential u* for mica surfaces with L z = 2.4<j and 
different particle numbers TV = 500, 1000, 2000. The symbols 
are the simulation data and the lines are a guide to the eye. 



rors. We conclude that the transitions found in figs. 0J) 
and |Ui are of first order. 

So far we have concentrated on density-driven crystal- 
lization transitions at two particular surface charges corre- 
sponding to silica and mica walls, respectively. However, 
as we will demonstrate below, such transitions also oc- 
cur when varying the surface charge continuously at fixed 
chemical potential (corresponding to a bulk density of 
p£ = 0.61). Our results for ijj n plotted in fig.|BJi reveal two 
major effects of such an "experiment", the first concern- 
ing the number of layers observed at given 0s- Indeed, 
increasing |0s| from zero one first finds a rather abrupt 
change (as revealed by density profiles) from a two- to a 
one-layer system, reflecting the effective narrowing of the 
pore due to increasing Coulomb repulsion (see also fig. [1} . 
Further charging then yields a second transformation back 
to a two-layer system, consistent with the fact that the 
range of repulsion decreases with increasing \ips\ for sur- 
face charges beyond |V>s| ~ 100mV. The charge-induced 
changes of the number of layers are accompanied by sud- 
den changes of the order parameter values, as one might 
expect. More significantly, however, is the fact that the or- 
der parameters exhibit further sudden changes within the 
ranges of constant layer number, accompanied by sharp 
maxima of the corresponding susceptibilities (see fig. Eb) . 
Specifically, increasing |0>s| from zero one detects a first 
transition (within the two-layer regime) at \ips\ ~ WmV, 
where the fluid crystallizes into a state with square-like 
order. For slightly larger |0s|, the system melts again 
into a fluid state. Similarly, there is a fluid-to-hexagonal 
transition at |-0s| ~ 80m V and a fluid-to-squarc transi- 
tion at \ips\ ~ IQOmV accompanied by large bond order 
susceptibilities. Note, that the transitions between one- 
and two-layer systems are not symmetric, i. e., the two- 
layer system "jumps" for small |0>s| into an one-layered 
fluid phase, whereas for high |0?s| it "jumps" into a one- 
layered hexagonal structure. Taken altogether, the results 
in fig. [5k, confirm the idea that variation of the surface 
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Fig. 6: Influence of surface charges (surface potential 0)s) on 
the in-plane order for L z = 2.4cr. (a) Bond order parameters 
•08 (squares) and tfig (triangles) as function of tps- The num- 
ber of layers are indicated by the vertical dashed lines and 
the numbers, (b) Corresponding bond order susceptibilities x& 
(squares) and \q (triangles), (c) Phase diagram as function of 
ips and bulk density p* h showing squared phases (squares, 'S') 
with 08 > 0.6, hexagonal phases (triangles, 'H') with -06 > 0.6, 
and non-ordered phases (circles, 'N'). The dashed lines are a 
schematic representation of the phase boundaries. 



charge can induce freezing and melting transitions. 

Finally, we give in fig. [U an overview of the system's 
behavior at the exemplary wall separation L z = 2.4<r in- 
cluding results at various densities in the range 0.57 < 
Pb — 0.614. The resulting "phase diagram" in the V's-Pb 
plane displays non-ordered (N), squared (S) and hexag- 
onal (H) phases. The diagram also contains the results 
of figs, [fjh, and [5b for p£ = 0.61. The observed phase se- 
quences between non-ordered and ordered phases includ- 
ing reentrance of the H phase are a result of the vary- 
ing shape of ups with -0s (see fig. QJ1) combined with the 
value of L z considered. The bond order threshold for dis- 
tinguishing non-ordered and ordered phases is set to 0.6, 
i. e., 06 > 0.6 is considered as hexagonal and 08 > 0.6 as 
squared. This criterion is motivated by the fact that the 
freezing transitions indicated by the peaks in the bond 
order susceptibilities Xn in fig- Hb correspond to values 
n « 0.6. However, the choice of the threshold does not 
change the qualitative behavior of the phase diagram. We 
note that at the parameters investigated we did not find 
any other complex, such as buckled, structures which have 
been observed in hard sphere systems [3J5] and ground- 
state Coulomb systems [3D]- One reason for the absence 
of buckling in our system might be the rather small value 
of k and the finite temperature. Indeed, according to a re- 
cent ground-state study [9], complex structures typically 
occur at larger screening parameters. Furthermore, a fi- 
nite temperature MC study of Yukawa bilayers [TU] found 
different (staggered) phases without buckling. In our sys- 
tem the region between one- and two-layer systems (see 
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fig- lit): where buckling might emerge, is characterized by 
a fluid-like order with small bond order parameters. 

Conclusion. — In summary, our computer simulation 
results demonstrate that surface charges have a crucial im- 
pact on the lateral ordering of charged colloidal suspen- 
sions confined to nanoscopic films. The observed effects 
originate from the impact of the wall counterions on the 
fluid- wall interaction, whereas the details of the fluid-fluid 
interaction play only a minor role. Most spectacularly, 
we find that charged walls support crystalline, hcxagonal- 
and square-ordered phases (depending on the actual sur- 
face potential) under conditions where the corresponding 
suspension between uncharged walls (and its bulk coun- 
terpart) is still fluid-like. Changing continuously the wall 
charge at fixed wall separation and density we find both, 
layering transitions and first-order, in-plane freezing and 
melting transitions including reentrant behavior. The ex- 
emplary cases studied in the present work indicate that 
charges on confining substrates add a wealth of structure 
formation phenomena to the already rich and complex 
phase diagram of confined charged colloids or dusty plas- 
mas. 

Of course, future experimental work is needed to con- 
firm these theoretical predictions, one first indication be- 
ing the unusual shape of the oscillatory surface forces 
characterizing dense suspensions of silica nanoparticles be- 
tween (charged) silica surfaces [TB]. Moreover, it would be 
very interesting and also important from a material science 
point of view [H2] to extend the investigations towards the 
case of binary colloidal mixtures of differently sized and 
charged particles. Work in this direction is under way |31j . 

We thank the Deutsche Forschungsgcmcinschaft for fi- 
nancial support via grant KL1215/6. 
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